PHS_Covid Analysis

Load the library

Load the dataset

Plot1(a): Trend on people who tested positive (For all the data).

#Plot to visualise trend on positive cases.
plot_positive_case <- trend_hb_daily %>% 
  filter (hb_name == "Scotland") %>% 
  ggplot()+
  aes(x = date, y = daily_positive)+
  geom_line()+
  scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  ggtitle("People tested positive") +
  xlab("Date (Months)") +
  ylab("No of Positive Cases") +
  color_theme() 

ggplotly(plot_positive_case)

NA
NA

Plot1(b): Trend on people who tested positive(For Past 6 months).

plot_positive_case <- trend_hb_daily %>% 
  filter(date >="2021-04-01") %>% 
  filter (hb_name == "Scotland") %>% 
  ggplot()+
  aes(x = date, y = daily_positive)+
  geom_line()+
  #scale_x_date(breaks = "1 month")+
  scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  ggtitle("People tested Positive") +
  xlab("Date (Months)") +
  ylab("No of Positive Cases") +
  color_theme()

ggplotly(plot_positive_case)

Plot1(c): Total Positive Case by neighborhood (Based on NHS Health board).

trend_hb_daily_tot <- trend_hb_daily %>% 
  filter(hb_name != "Scotland") %>% 
  group_by(hb) %>% 
  summarise(total_count = sum(daily_positive))
  

map_data <- zones_nhs_hb %>% 
    inner_join(trend_hb_daily_tot, by = c("code" = "hb"))
  
map_palette <- colorNumeric("viridis", domain = range(map_data$total_count))

map_data <- map_data %>% 
  arrange(total_count)

 map_data <- map_data %>% 
    mutate(colour = map_palette(total_count))
 
  map_data %>%
    leaflet() %>%
    addPolygons(
      popup = ~ str_c(name, "<br>",  "Positive Case = ", total_count, sep = ""),
      color = ~colour
    ) %>%
    addLegend(
      position = "bottomright",
      colors = ~colour,
      labels = ~name
    )

NA

Plot1(c): Total Positive Case by neighborhood (Based on Local Authority).

trend_la_daily_tot <- trend_la_daily %>% 
  group_by(ca) %>% 
  summarise(total_count = sum(daily_positive))

map_data <- zones_la %>% 
    inner_join(trend_la_daily_tot, by = c("code" = "ca"))
  
map_palette <- colorNumeric("viridis", domain = range(map_data$total_count))

map_data <- map_data %>% 
  arrange(total_count)

 map_data <- map_data %>% 
    mutate(colour = map_palette(total_count))
 
map_data %>%
    leaflet() %>%
    addPolygons(
      popup = ~ str_c(name, "<br>",  "Positive Case = ", total_count, sep = ""),
      color = ~colour
    ) %>%
    addLegend(
      position = "bottomright",
      colors = ~colour,
      labels = ~name
    )

EXTRA 1: Forecast on Positive cases:

# The data is chosen from Jun as it is the beginning of wave 1 
trend <- trend_hb_daily %>% 
  filter (hb_name == "Scotland") %>% 
  filter(date >="2021-06-01") %>% 
  select(date, daily_positive)

# Data is converted into timeseries for forecasting
trend <- as_tsibble(trend, index = date)

# Different forecasting methods are used 
fit <- trend %>%
  model(
    snaive = SNAIVE(daily_positive),
    mean_model = MEAN(daily_positive),
    arima = ARIMA(daily_positive)
  )

# Forecast for 2 weeks (14 days)
forecast_14days <- fit %>%
  forecast(h = 14)
forecast_14days

forecast_14days %>%
 #filter(.model == "snaive") %>%
 autoplot(trend, level = NULL) +
  ggtitle("Forecasts for Positive cases for 2 weeks") +
  xlab("Year") +
  ylab("No of Positive Cases") +
  guides(colour = guide_legend(title = "Forecast"))+
  scale_x_date(breaks = "1 month", date_labels = "%B" )+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

EXTRA: Forecast for one month:

forecast_1month <- fit %>%
  fabletools::forecast(h = "1 month")

forecast_1month %>%
 #filter(.model == "snaive") %>%
 autoplot(trend, level = NULL) +
  ggtitle("Forecasts for Positive cases for one month") +
  xlab("Year") +
  ylab("No of Positive Cases") +
  guides(colour = guide_legend(title = "Forecast"))+
  scale_x_date(breaks = "1 month", date_labels = "%B" )+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Check if the forecast data fits our existing data

# Set the training data from Jan to Aug
train <- trend %>%
  filter_index("2021-01-01" ~ "2021-08-31")

# run the model on the training set 
fit_train <- train %>%
  model(
    mean_model = MEAN(daily_positive),
    arima = ARIMA(daily_positive),
    snaive = SNAIVE(daily_positive)
  )

forecast_test <- fit_train %>% 
  forecast(h = 30)

forecast_test %>%
  autoplot(train, level = NULL) +
    autolayer(filter_index(trend, "2021-09-01" ~.), color = "black") +
    ggtitle("Forecasts for positive cases") +
    xlab("Date") + ylab("Daily Positive") +
    guides(colour=guide_legend(title="Forecast"))

Plot variable not specified, automatically selected `.vars = daily_positive`

Check for the accuracy of the model

The minimum RMSE is chosen.

accuracy_model <- accuracy(forecast_test, trend)

accuracy_model %>% 
  select(-.type) %>%
  arrange(RMSE)

2 Analyse the trend on Hospitalisations:

plot_hosp <- trend_hb_daily %>% 
  filter(date >="2021-04-01") %>% 
  filter (hb_name == "Scotland") %>% 
  filter(!(is.na(hospital_admissions))) %>% 
  ggplot()+
  aes(x = date, y = hospital_admissions)+
  geom_col()+
  scale_x_date(breaks = "1 week", date_labels = "%d - %b" )+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  ggtitle("Trend on Hospitalisations") +
  xlab("Date (months)") +
  ylab("No of patients admitted") +
  color_theme()

ggplotly(plot_hosp)

NA

2 (b) Forecast on Hospitalisation:**

Using 6 months of data ( Two - Wave)

trend_hosp <- trend_hb_daily %>% 
  filter (hb_name == "Scotland") %>% 
  filter(date >="2021-06-01") %>% 
  filter(!(is.na(hospital_admissions))) %>% 
  select(date, hospital_admissions)

trend_hosp <- as_tsibble(trend_hosp, index = date)

fit <- trend_hosp %>%
  model(
    snaive = SNAIVE(hospital_admissions),
    mean_model = MEAN(hospital_admissions),
    arima = ARIMA(hospital_admissions),
    ets = ETS(log(hospital_admissions) ~ error("M") + trend("Ad") + season("A"))
  )

forecast_14days <- fit %>%
  forecast(h = 14)

forecast_14days %>%
# filter(.model == "snaive") %>%
 autoplot(trend_hosp, level = NULL) +
  ggtitle("Forecasts for Hospital admissions cases for 2 weeks") +
  xlab("Year") +
  ylab("No of patients admitted") +
  guides(colour = guide_legend(title = "Forecast"))+
  scale_x_date(breaks = "1 month", date_labels = "%B" )+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

forecast_1month <- fit %>%
  fabletools::forecast(h = "1 month")

forecast_1month %>%
 filter(.model %in% c("snaive","arima")) %>%
 autoplot(trend_hosp, level = NULL) +
  ggtitle("Forecasts for Hospitalisation for one month") +
  xlab("Year") +
  ylab("No of Positive Cases") +
  guides(colour = guide_legend(title = "Forecast"))+
  scale_x_date(breaks = "1 month", date_labels = "%B" )+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

#  set the training data
train <- trend_hosp %>%
  filter_index("2021-06-01" ~ "2021-08-31")

# run the model on the training set 
fit_train <- train %>%
  model(
    mean_model = MEAN(hospital_admissions),
    arima = ARIMA(hospital_admissions),
    snaive = SNAIVE(hospital_admissions),
    ets = ETS(log(hospital_admissions) ~ error("M") + trend("Ad") + season("A"))
  )

forecast_test <- fit_train %>% 
  forecast(h = 35)

forecast_test %>%
  filter(.model %in% c("arima","snaive")) %>%
  autoplot(train ) +
    autolayer(filter_index(trend_hosp, "2021-06-01" ~.), color = "black") +
    ggtitle("Forecasts for Hospitalisations") +
    facet_wrap(~.model)+
    xlab("Date") + 
    ylab("No of Patients admitted") +
    guides(colour=guide_legend(title="Forecast"))

Plot variable not specified, automatically selected `.vars = hospital_admissions`

accuracy_model <- accuracy(forecast_test, trend_hosp)

Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing. 
1 observation is missing at 2021-10-05

accuracy_model %>% 
  select(-.type) %>%
  arrange(RMSE)

Plot3(a): Trend on people who tested positive (For all the data).

#Plot to visualize trend on positive cases.
plot_vaccine <- daily_vacc_hb_plot %>% 
  ggplot()+
  aes(x = date, y = number_vaccinated)+
  geom_line(aes(color = dose))+
  facet_wrap(~dose)+
  scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  ggtitle("Trend on Vaccination") +
  xlab("Year") +
  ylab("No of Positive Cases") +
  color_theme()

ggplotly(plot_vaccine)